Stoffbilanzen (N)
$\dot n_i = \dot n_{i, 0} + \mathcal{R} \dot n_{i, v} + \sum_{j}{\nu_ij \xi_j}$
Energiebilanz (1)
$\begin{array}{lll} 0 &= \dot Q & + \sum\limits_i(\dot n_i (\Delta H_i^{\circ}(T)-\Delta H_{i,0}^\circ))_{ein}- \sum\limits_i(\dot n_i (\Delta H_i^{\circ}(T)-\Delta H_{i,0}^\circ))_{aus} + \sum\limits_{j}{\xi_j (-\Delta Hr_j(T))}\\ &= 0 &+ \sum\limits_i((\dot n_{i,0}+\mathcal{R}\cdot\dot n_{i,v}) \cdot(\Delta H_i^{\circ}(T_{ein})-\Delta H_{i,0}^\circ))\\ &&- \sum\limits_i(\dot n_{i,2}\cdot(\Delta H_i^{\circ}(T_2)-\Delta H_{i,0}^\circ)) \\ &&+\sum\limits_{j}{\xi_j \cdot (-\Delta Hr_j(T))}\\ \end{array}\\ $
Gleichgewichtskonstanten (r)
$\begin{array}{ll} K_j(T) &= exp \left(-\frac{\Delta H_0^\circ}{R T} + \frac{(\Delta H_0^\circ -\Delta G_0^\circ)}{R T_0} - \frac{1}{T}\int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} + \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}\right) \\ &= \prod_i (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum_i \nu_i}(n)^{-\sum_i \nu_i}\end{array}$
$p^0 = 1 bar$
Idealer Gas-Ansatz, $K_{\phi^{eq}}=1$
Stoffbilanzen (N)
$\dot F z_i = \dot V y_i + \dot L x_i $
Isotherme Verdampfungs-Gleichgewichtsbedingungen (N)
$(\hat{\phi_i}^V P) y_i^V=(\hat{\phi_i}^L P) x_i^L \hspace{2cm} K_i\equiv{{y_i}\over{x_i}}={{\hat{\phi_i}^L}\over{\hat{\phi_i}^V}}$
Zusatzgleichung (1)
$0 = \sum\limits_i y_i - \sum\limits_i x_i$
Verringerung des Systems zu einer Funktion einer Veränderlichen: $\psi\equiv(V/F)$
$\begin{array}{ll} z_i &= \psi y_i + (1-\psi) x_i \\ &= \psi x_i K_i + (1-\psi) x_i \\ x_i &={{z_i}\over{1+\psi(K_i-1)}} \hspace{3cm} y_i &={{z_i K_i}\over{1+\psi(K_i-1)}}\\ \Rightarrow 0 &= \sum\limits_i y_i - \sum\limits_i x_i = \sum\limits_i {{z_i(K_i-1)}\over{1+\psi(K_i-1)}}\\ \end{array}$
Stoffbilanzen (N)
$\dot V y_i = \dot n_{i,v} = \dot n_{i,r} + \dot n_{i,p}$
Splitter-Bedingungen (N)
$\dot n_{i,r} = \mathcal{R} \dot n_{i,v}$
$\Rightarrow \dot n_{i,p} = (1-\mathcal{R}) \dot n_{i,v}$
Eingangs Daten: Zusammensetzungen der Gasphase, der flüssigen Phase, Temperatur und Druck $\{x_i\}, \{y_i\}, T, P$
Datenausgabe: Realgasfaktor beider Phasen $Z^L, Z^V$ Fugazitätskoeffizienten der Komponenten in beiden Phasen $\{\hat{\phi_i^L}\}, \{\hat{\phi_i^V}\}$, Verteilungskoeffizienten $\{K_i\}$.
Vorgehensweise:
Phasen-unabhängige Terme berechnen
$\begin{array}{|c|c|} \hline a_i(T) = \Psi{{ \alpha(Tr_i, \omega_i) R^2 Tc_i^2}\over{Pc_i}}& b_i =\Omega{{R Tc_i}\over{Pc_i}}\\ \hline \beta_i = {{b_i P}\over{R T}}& q_i = {{a_i(T)}\over{b_i R T}} \\ \hline \end{array}$
Phasen-abhängige Mischregeln für jede Phase berechnen (L, V)
$\begin{array}{|c|c|} \hline a = \sum\limits_i\sum\limits_j x_i x_j a_{ij} & b=\sum\limits_i x_i b_i\\ a_{ij} = (a_i a_j)^{1/2} & \\ \hline \beta^p \equiv {{b^p P}\over{R T}} \hspace{10mm} \textit{p: (L, V)} & q^p \equiv {{a^p}\over{b^p R T}} \hspace{10mm} \textit{p: (L, V)} \\ \hline \end{array}$
Partielle molare Größen in jeder Phase berechnen
$\require{cancel} \begin{array}{|ccll|} \hline \bar{a_i} &\equiv\left[{{\partial (n a)}\over{\partial n_i}} \right]_{T,n_j} & =\left[{{\partial }\over{\partial n_i}}(n \sum\limits_i\sum\limits_j x_i x_j a_{ij}) \right]_{T,n_j} &=-a + 2\sum\limits_{j\neq i} x_j a_{ij} + 2 x_i a_i\\ \bar{b_i} &\equiv\left[{{\partial (n b)}\over{\partial n_i}} \right]_{T,n_j} &= \left[{{\partial (n_i b_i)}\over{\partial n_i}} \right]_{T,n_j} + \cancelto{0}{\sum\limits_{j\neq i}\left[{{\partial (n_j b_j)}\over{\partial n_i}} \right]_{T,n_j}} &=b_i \\ \bar{q_i} &\equiv \left[{{\partial (n q)}\over{\partial n_i}} \right]_{T,n_j} &= q \left(1+{{\bar{a_i}}\over{a}}-{{\bar{b_i}}\over{b}}\right) &= q \left(1+{{\bar{a_i}}\over{a}}-{{b_i}\over{b}}\right)\\ \hline \end{array}$
Nach den Realgasfaktor $Z^L$ bzw. $Z^V$ jeder Phase die Kubische Zustandsgleichung auflösen:
Für den Gas ist die "V Form" am geeignetesten, mit erster Schätzung $Z=1,0$. Wenn es sich um einem tatsächlichen Gleichgewichtszustand handelt, soll die Lösung Gas-ähnlich (höhere Nullstelle) sein.
$\begin{array}{|ll|} \hline L \hspace{2cm} Z^L &= \beta^L+(Z^L+\epsilon \beta^L)(Z^L+\sigma \beta^V)\left({{1+\beta^L-Z^L}\over{q^L \beta^L}}\right) \\ V \hspace{2cm} Z^V &= 1 + \beta^V - q^V \beta^V {{Z^V-\beta^V}\over{(Z^V+\epsilon \beta^V)(Z^V+\sigma \beta^V)}}\\ \hline \end{array}$
Fugazitätskoeffizienten und Verteilungskoeffizienten Mit den Lösungen des Realgasfaktors berechnen sich die Fugazitätskoeffizienten.
$\begin{array}{|ll|} \hline ln(\hat{\phi_i})&={{b_i}\over{b}}(Z-1)-ln(Z-\beta)-\bar{q_i} I\\ I &={{1}\over{\sigma-\epsilon}}ln\left({{Z+\sigma \beta}\over{Z+\epsilon \beta}} \right)\\ K_i &={{\hat{\phi_i}^L}\over{\hat{\phi_i}^V}}\\ \hline \end{array}$
In der Dampfphase mit den oben bestimmten Parametern, Residualeigenschaften berechnen
$\begin{array}{|ll|} \hline \frac{G^R}{R T} &= Z-1-ln(Z-\beta)-q I\\ \frac{H^R}{R T} &= Z-1+\left[ \frac{d\thinspace ln \thinspace \alpha(T_r)}{ d \thinspace ln \thinspace T_r} -1\right]q I\\ \frac{S^R}{R} &= ln(Z-\beta)+ \frac{d\thinspace ln \thinspace \alpha(T_r)}{ d \thinspace ln \thinspace T_r} q I\\ \hline \end{array}$
Da $M^{R}\equiv M - M^{Id}$, anhand der aus tabellierten Korrelationen zu ermittelnden Idealgaseigenschaften, die Realeigenschaften berechnen.
$\begin{array}{ll} H^{Id} &= H_0^{Id} + \int\limits_{T_0}^{T}{Cp^{Id}dT}\\ S^{Id} &= S_0^{Id} + \int\limits_{T_0}^{T}{\frac{Cp^{Id}}{T}dT}- R \thinspace ln(\frac{p}{p^0})\\ G^{Id} &= H^{Id} - T S^{Id}\\ \end{array}$
Ausdrücke einschließlich Idealgasanteil und Residualeigenschaften:
$\begin{array}{|ll|} \hline H &= H_0^{Id} + R \cdot \int\limits_{T_0}^{T}{\frac{Cp^{Id}}{R}dT} + H^R\\ S &= S_0^{Id} + R \cdot \int\limits_{T_0}^{T}{\frac{Cp^{Id}}{R T}dT}- R \thinspace ln(\frac{p}{p^0})+S^R\\ G &= H - T S\\ \hline \end{array}$
$\begin{array}{ll} exp \left(- \frac{\Delta G_i}{R T} \right) &= K_p K_{\phi^{eq}} = K_x \prod\limits_{i} \left( \frac{p}{p^0}\right)^{\nu_i} K_{\phi^{eq}} \\ &=\prod\limits_{i} (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum\limits_{i} \nu_i}(n)^{-\sum\limits_{i} \nu_i} K_{\phi^{eq}}\end{array}$
$\Delta H^\circ = \Delta H_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT}$
$\Delta S^\circ = \Delta S_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$
$\Delta G^\circ = \Delta H^\circ - T \Delta S^\circ = \Delta H_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} - T \Delta S_0^\circ - R T \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$
$\Delta S_0^\circ = \frac{\Delta H_0^\circ - \Delta G_0^\circ}{T_0}$
$\Delta G^\circ = \Delta H_0^\circ - \frac{T}{T_0}(\Delta H_0^\circ -\Delta G_0^\circ) + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} - R T \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$
$\begin{array}{ll} K_{(T)} &= exp \left(-\frac{\Delta H_0^\circ}{R T} + \frac{(\Delta H_0^\circ -\Delta G_0^\circ)}{R T_0} - \frac{1}{T}\int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} + \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}\right) \\ &= \prod_i (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum_i \nu_i}(n)^{-\sum_i \nu_i}\end{array}$
Somit läßt sich K(T) bestimmen, insofern man über einen Ausdruck für $Cp_i(T)$ verfügt. Bei geringer Veränderung der Wärmekapazität Cp im Temperatur-Bereich kann man auch einen bestimmten Mittelwert als ~konstant einsetzen.
In [4]:
from scipy import optimize
import numpy as np
p = 50. # bar
temp = 273.15 + 220. # K
t_flash = 273.16 + 60 # K
t0_ref = 298.15 # K
r = 8.314 # J/(mol K)
rvg = 0.8 # Rückvermischungsgrad
namen = ['CO', 'H2', 'CO2', 'H2O', 'CH3OH', 'N2']
n0co = 750. # kmol/h
n0h2 = 5625. # kmol/h
n0co2 = 750. # kmol/h
n0h2o = 375. # kmol/h
n0ch3oh = 0. # kmol/h
n0n2 = 500. # kmol/h
ne = np.array([n0co, n0h2, n0co2, n0h2o, n0ch3oh, n0n2])
nuij = np.array([[-1, -2, 0, 0, +1, 0] ,
[0, -3, -1, +1, +1, 0],
[-1, +1, +1, -1, 0, 0]]).T
h_298 = np.array(
[-110.541, 0., -393.505, -241.826,-201.167, 0.]) * 1000 # J/mol
g_298 = np.array(
[-169.474, -38.962, -457.240, -298.164, -272.667, -57.128]) * 1000 # J/mol
# Berechne delta Cp(T) mit Temperaturfunktionen für ideale Gase (SVN).
# Koeffizienten für Cp(T)/R = A + B*T + C*T^2 + D*T^-2, T[=]K
# Nach rechts hin: A, B, C, D
# Nach unten hin: CO, H2, CO2, H2O, CH3OH, N2
cp_coefs = np.array([
[
y.replace(',', '.') for y in x.split('\t')
] for x in """
3,3760E+00 5,5700E-04 0,0000E+00 -3,1000E+03
3,2490E+00 4,2200E-04 0,0000E+00 8,3000E+03
5,4570E+00 1,0450E-03 0,0000E+00 -1,1570E+05
3,4700E+00 1,4500E-03 0,0000E+00 1,2100E+04
2,2110E+00 1,2216E-02 -3,4500E-06 0,0000E+00
3,2800E+00 5,9300e-04 0,0000E+00 4,0000e+03
""".split('\n') if len(x)>0], dtype=float)
def cp(t):
return r * (
cp_coefs[:,0] +
cp_coefs[:,1] * t +
cp_coefs[:,2] * t**2 +
cp_coefs[:,3] * t**-2
) # J/(mol K)
# Berechne H(T), G(T) und K(T) mit Cp(T)
def h(t):
return (
h_298 +
r * cp_coefs[:,0]*(t-t0_ref) +
r * cp_coefs[:,1]/2.*(t**2-t0_ref**2) +
r * cp_coefs[:,2]/3.*(t**3-t0_ref**3) -
r * cp_coefs[:,3]*(1/t-1/t0_ref)
) # J/mol
def g(t, h_t):
return (
h_t - t/t0_ref*(h_298 - g_298) -
r * cp_coefs[:,0]*t*np.log(t/t0_ref) -
r * cp_coefs[:,1]*t**2*(1-t0_ref/t) -
r * cp_coefs[:,2]/2.*t**3*(1-(t0_ref/t)**2) +
r * cp_coefs[:,3]/2.*1/t*(1-(t/t0_ref)**2)
) # J/mol
def k(t, g_t):
delta_g_t = nuij.T.dot(g_t)
return np.exp(-delta_g_t/(r * t))
delta_gr_298 = nuij.T.dot(g_298)
delta_hr_298 = nuij.T.dot(h_298)
cp_493 = cp(493.15) # J/(mol K)
h_493 = h(493.15) # J/mol
g_493 = g(493.15, h_493) # J/mol
k_493 = k(493.15, g_493) # []
for i, f in enumerate(delta_hr_298):
print('Delta H_' + str(i+1) + '(298.15K)=' + str(f/1000.) + 'kJ/mol')
print('\n')
for i, f in enumerate(k_493):
print('K' + str(i+1) + '(493K)=' + str(f))
print('\n')
n0 = np.array([n0co, n0h2, n0ch3oh])
def fun(x_vec):
nco = x_vec[0]
nh2 = x_vec[1]
nch3oh = x_vec[2]
xi1 = x_vec[3]
t = x_vec[4]
n = np.array([nco, nh2, nch3oh])
cp_t = cp(t)
h_t = h(t)
g_t = g(t, h_t)
k_t = k(t, g_t)
h_ein = h_t[[1, 2, -2]]
cp_ein = cp_t[[1, 2, -2]]
cp_t = cp_t[[1, 2, -2]]
h_t = h_t[[1, 2, -2]]
g_t = g_t[[1, 2, -2]]
delta_h_t = nuij[[1, 2, -2]].T.dot(h_t) # J/mol
f1 = -nco + n0co - xi1
f2 = -nh2 + n0h2 -2*xi1
f3 = -nch3oh + n0ch3oh +xi1
f4 = -k_t[0] * (nco * nh2**2) + \
nch3oh * (p/1.)**-2 * (nco + nh2 + nch3oh)**-(-2)
f5 = np.sum(
np.multiply(n0, cp_ein)*temp -
np.multiply(n, cp_t)*t
) + xi1 * (-delta_h_t[0])
return [f1, f2, f3, f4, f5]
x0 = np.append(n0, [0., temp])
sol = optimize.root(fun, x0)
f_final = - sol.x[:3].reshape([3,1]) + ne[[0,1,4]].reshape([3,1]) + nuij[:,0][[0,1,4]].reshape([3,1])*sol.x[-2]
print(sol)
print('\n\n')
print('Zustand der Optimisierungs-Funktionen\n')
print(f_final)
print('\n\n')
print('T_ein=493.15K, p=50 bar, in adiabatischem Reaktor')
print('Lösung für nur einzige Reaktion (ohne CO2):\n')
for i, f in enumerate(sol.x[:2]):
print('n_' + namen[i] + '= ' + str(f) + ' kmol/h')
print('n_' + namen[4] + '= ' + str(sol.x[2]) + ' kmol/h')
print('T= ' + str(sol.x[-1]) + ' K')
n0 = np.array([n0co, n0h2, n0co2, n0h2o, sol.x[2], n0n2])
#n0 = ne
# Lösung des einfacheren Falls in schwierigerem Fall einwenden.
def fun(x_vec):
nco = x_vec[0]
nh2 = x_vec[1]
nco2 = x_vec[2]
nh2o = x_vec[3]
nch3oh = x_vec[4]
nn2 = x_vec[5]
xi1 = x_vec[6]
xi2 = x_vec[7]
xi3 = x_vec[8]
t = x_vec[9]
n = np.array([nco, nh2, nco2, nh2o, nch3oh, nn2])
xi = np.array([xi1, xi2, xi3])
h_ein = h_493
cp_ein = cp_493
cp_t = cp(t)
h_t = h(t)
g_t = g(t, h_t)
k_t = k(t, g_t)
delta_h_t = nuij.T.dot(h_t) # J/mol
f1 = -nco + n0co - xi1 +0 -xi3
f2 = -nh2 + n0h2 -2*xi1 -3*xi2 +xi3
f3 = -nco2 + n0co2 +0 -xi2 +xi3
f4 = -nh2o + n0h2o +0 +xi2 -xi3
f5 = -nch3oh + n0ch3oh +xi1 +xi2 -0
f6 = -nn2 + n0n2 + 0
f7 = -k_t[0] * (nco * nh2**2) + \
nch3oh * (p/1.)**-2 * (nco + nh2 + nco2 + nh2o + nch3oh + nn2)**-(-2)
f8 = -k_t[1] * (nco2 * nh2**3) + \
nch3oh * nh2o * (p/1.)**-2 * (nco + nh2 + nco2 + nh2o + nch3oh + nn2)**-(-2)
f9 = -k_t[2] * (nco * nh2o) + \
nco2 * nh2 * (p/1.)**0 * (nco + nh2 + nco2 + nh2o + nch3oh + nn2)**-0
f10 = np.sum(
np.multiply(n0, (h_ein-h_298)) -
np.multiply(n, (h_t-h_298))) + np.dot(xi, -delta_h_t)
return [f1, f2, f3, f4, f5, f6, f7, f8, f9, f10]
x0 = np.append(n0, [0., 0., 0., sol.x[-1]])
sol = optimize.root(fun, x0)
print('\n\n')
print('success: ' + str(sol.success))
f_final = - sol.x[:6].reshape([6,1]) + ne.reshape([6,1]) + nuij.dot(sol.x[6:-1].reshape([3,1]))
print('\n\n')
print('T_ein=493.15K, p=50 bar, in adiabatischem Reaktor.')
print('Lösung für alle drei Reaktionen, mit CO2:\n')
for i, f in enumerate(sol.x[:6]):
print('n_' + namen[i] + '= ' + str(f) + ' kmol/h')
print('\n')
for i, f in enumerate(sol.x[6:-1]):
print('xi_' + str(i) + '= ' + str(f) + ' kmol/h')
print('\n')
print('T=' + str(sol.x[-1]) + ' K, oder...')
print('T=' + str(sol.x[-1]-273.15) + ' °C')
print('\n')
print('0 = Q + Sum(Delta H)_ein - Sum(Delta H)_aus')
bilanz = np.sum(
np.multiply(n0, (h_493-h_298)) -
np.multiply(sol.x[:6], (h(sol.x[-1])-h_298))
) + np.dot(sol.x[6:-1], -nuij.T.dot(h(sol.x[-1])))
annaeherung = np.sum(
np.multiply(n0, cp_493)*493.15 -
np.multiply(sol.x[:6], cp(sol.x[-1]))*sol.x[-1]
) + np.dot(sol.x[6:-1], -nuij.T.dot(h(sol.x[-1])))
print('-Q = (n.(H_t-H_298))_ein -(n.(H_t-H_298))_aus + Sum(xi_j * (-Delta Hr_j)) = ' +
str(bilanz) + 'J/h')
print('\n\n')
print('Zustand der Optimisierungs-Funktionen\n')
print(f_final)
print('\n\n')
print('Umsatz(CO): ' +
'{:.4g}'.format((ne[0]-sol.x[0])/ne[0]))
print('Umsatz(CO2): ' +
'{:.4g}'.format((ne[2]-sol.x[2])/ne[2]))
print('Niedriger Umsatz!')
print('Ausbeute (CH3OH/CO): ' +
'{:.4g}'.format((sol.x[4]-ne[4])/(ne[0]-sol.x[0])))
print('Ausbeute (CH3OH/CO2): ' +
'{:.4g}'.format((sol.x[4]-ne[4])/(ne[2]-sol.x[2])))
In [8]:
import z_l_v
z_l_v.beispiel_pat_ue_03_vollstaendig(0.65, print_output=True)
Out[8]: